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source code implicit finite element solver 



ra te/temp/press dependents von mises isotropic plasticity 
umat for abaqus S-S- nonlinear strain hardening. 
Hd/3d problems with the exception of plane stress 
by omar a hasan last modified DS-Oe-Tb 

user must specify differential hardening data in umat 
and dimension hardening table appropriately 
must have atleast two sets of points in table 
subroutine umatCstressTStateviddsddeTSseTspdiScd-i 
1 rplnddsddtidrpldeidrpldti 

B stranTdstranTtimeidtime-itempTdtempTpredefTdpredTCninamGT 
3 ndi -.nshr-.ntenSTnstatvnpropSinpropSTCoordSidrot-ipnewdt-i 
M celentidfgrdD-.dfgrdl'.noelinptTlayerikspt-.kstepTkinc) 



include ' aba__par am . i nc ' 
character*a cmname 

dimension strossCntons) - statevCnstatv) -i 
1 ddsddeCntensTntens)Tddsddt(ntens)idrplde(ntens)T 
a strdn(ntens)-idstran(ntens)'.time(5)Tpredef(l)ipred(l)T 
3 props (npr ops) iCoords(3)-.drot(3-i3)->dfgrdQ<3-.3) d f grdl ( 3 -» 3) 

dimension flowCb) 

parameter ( zero = a . dO Tone = l * dO t two = e . dO three = 3 . daTsix = t, . dQ^ 
1 newton = ba toler = l . Dd-S-i twbth = D . bbbbbbbbbfeifcidD ) 



cannot be used for plane stress 



props(l) - e (Pa) (temperature dependent) 
props(2) - nu 

props(3) - rate sensitivity (temperature dependent) 
props(H) - intrinsic flow rate (temperature dependent) 
props(5) - pressure sensitivity 

calls uhard for curve of intrinsic strength vs. plastic strain 



material properties 
e.mod = propsCl) 
enu=props ( 2 ) 

ebulk3=efflod/(one-two*enu) 

eg2=emod/(one+enu) 

eg = eg2/ 1 wo 

eg3=three*eg 

elam=(ebuik3-eg2)/three 

rlp2m=elam+eg2/three 

rates f =pr ops < 3 ) 

rrates=one/ratesf 

dtebsa=dtime*props(M) 
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ps f = props ( S ) 

esi=dstran(l)**^+d5tran(2)**2+dstran(3)**2 
do kl = ndi + l -tntens 

esi = esi + tuo*(dstran(kl)/tiJo)**B 
end do 

esi=sqrt(twbth*esi) 

s_rate = max (l-d-lQ-iesi/dtime) 

c 

c elastic stiffness 

call aset(ddsdde-.zero-tntens*ntens) 
do kl = l-.ndi 
do kE^lnndi 

ddsddeCkankD^elam 
end do 

ddsdde(klnkl)=eg2+elam 
Qnd do 

do kl=ndi+l T ntsns 

ddsdds (kl^kl) =eg 
end do 

c 

c recover equivalent plastic strain a equivalent stress 
c and hydrostatic stress at start of step 

eqplas=statev (1) 

qold=statev ( H) 

hydr_o = (stress(l)+5tre5sC2)-»-stress(3))/three 

c 

c calculate predictor stress 
do kl=lnntens 
do k2 = l T n tens 

stress (ka)=str ess Ck2) + ddsdde (kE ■.kl)*dstr an (kl) 

end do 
end do 

c 

c calculate equivalent von mises stress 

smises=Cstress<l)-stress(E)) ( stress ( E ) -stress ( 3 ) )**E 
I +(stress(3)-stress(l) )**E 

do kl = ndi + l-» ntens 

smises=smises+six*stress (kl) **E 
end do 

smises=sqrt(smises/two) 

c 

c get differential hardening from the specified hardening curve 
call uhard(syielO-.hardieqplas) 

c 

c determine if actively yielding 
if (timeCl) .gt-O-dO) then 

c 

c separate the hydrostatic from the deviatoric stress 

c calculate the flow direction 



Cantor CoL. ..i LLP 55 Griffin Road South, Bloomfield, CT 06002 
(860) 2S6-2929 



source code implicic finite element: solver 



shyclro=(stress(l)+stress(5)+stress(3))/three 

do l<l = lTndi 

flow(kl)=(strGSS(kl)-shydro)/5mises 

end do 

do kl = nd i + 1 intens 

flow(kl)=stress(kl)/smises 
end do 

solve for equivalent von mises stress 

and equivalent plastic strain increment using newton iterati 
sy ield=sy ielO 

use this to minimize iterations during elastic deformation C 

dGqpl=dtebsD*exp((smises-syiGld)*ratesf) 
use this to minimize iterations during plastic deformation ( 

deqpl =es i 

do kewton=l ineuton 

deqpl=fflax(deqplnl-d-50) 

qhs=smises-Gg3*deqpl-syiGld-rrates*dlogCdeqpl/dtebsQ) 
rhs=qhs+psf*shydro 

deqpl=deqpl+deqpl*rhs/Cdeqpl*(eg3+hard)+rrates) 
call uhardCsyieldihardneqplas+deqpl) 
if (abs(rhs) .lt.toler*bD.dD) goto 10 
end do 

write(7i2) newton 

format (// iBDx-. ' ***warning - plasticity algorithm did not 

'converge after 'ii3-i' iterations') 
write(7-.*)dstran(l)-,dstran(2)'.dstran(3)-.dstran(M) 
u)rite(7T*)dstran(S)-.dstran(b) -.esiTsmises-iStatevCl) 
write(7'.*)statev(2)-.statev(3)-.statevCM)-.statev<5) 
write(7T*)qhSTdeqplTrhSTShydroi5tress(l)nStress(E) 
writGC7-i*)stress(3)TStress(M)TStres5(b)TStressCb) 
continue 

the new equivalent deviatoric stress (q) is 
q=syield+rrates*dlog (deqpl/dtebsD) -psf*shydro 

update stressT elastic and plastic strains and 
equivalent plastic strain 

do kl = l T nd i 

stress(kl)=flow(kl)*q+shydro 

end do 

do kl=ndi+lTntGns 

stress(kl)=flow(kl)*q 
end do 

eqplas=eqplas+deqpl 
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calculate plastic dissipation 
spd=deqpl*(qold+q) /two 

formulate the jacobian (material tangent) 
first calculate effective moduli 
e f f g=eg*q/sm i ses 
ef f gE=two*ef f g 
effg3=three/tuo*effgE 
effiam=<ebulk3-effg5)/three 
hardl=hard+rrates/deqpl 
effhrd=eg3*hardl/(eg3+hardl)-effg3 
cee=-ebulk3*psf*eg*deqpl/smises 
do kl = l-indi 
do k2 = linndi 

ddsdde(ka'.kl)=efflam + cee*flow(k5) 
end do 

ddsdde(kl'.kl) = effgE + efflam + cee*flow(kl) 
end do 

do kl=ndi+lTntens 

ddsdde(klTkl)=ef f g 
end do 

do kl=lintens 
do k2 = l-tntens 

ddsdde(k2-.kli)=ddsdde(k2-.kl)+effhrd*flow(k2)*flow(kl) 



end do 
endi f 

store state variables in arrs'^ 

equiv s tra i n mises s tress p 1 as t i c strain rate^elastic strain 
rate and iterations to convergence 

statev Cli)«eqplas 

statev (E) =q 

statev (3)=deqpl/dtime 

st atev ( M ) =es i/d tima 

statev (S) =kewton 

return 
end 

subroutine uhardCsyieldihardTeqplas) 
include ' aba_param . inc ' 

table must be dimensioned correctly below: 
dimension table (En?) 
par ame ter ( zer 0 = 0 - dO ) 
nbv 313 hardening table 
n va lue = 7 

this is room temp data 
tabled^D^O-aOdO 



end do 
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tableCEil) =0-0 
tablediE) ==-5 .a^SdO 
table(E-.E)=0.1Sl 
table(l-.3) =-3 .OMdO 
table(5.3)=D.337 
table(l-.4)=M.75Lda 
table(E.M) SME 
table (l.S) =m . MldD . 
tableCETS)=0.73ti 
tableCl-ib)=4a.l4bdQ 
tabLeCE-.b)=l-aT3 
tabled. 7)=E7a4.MdO 
table(ai7)=17.aflb 

do kl=l T nva 1 u e-1 

Qqpll = tablQC2-.k3j + l) 
ifCeqplas-lt-eqpll) then 
eqplO=table(Enkl) 
current yield stress and hardening 

deqpl=eqpll-eqplO 
syielO=table (likl) 
syiell = table (l-.kl + l) 
dsyiel=syiell-syielO 
hard=dsyiel/deqpl 

5yield=syielD+(eqplas-eqplD)*hard 

goto 10 
endi f 
end do 
continue 

return 
end 
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c vectorized user material subroutine for shell and plane 

c stress elements (abaqus5-S) 

c rate/temp dependent isotropic plasticity with linear 

c elasticity! strain softening/hardening a press- depnd. 

c yield 

c by omar a hasan ( has anSlcr d . ge ■ com ) 

c last modified □S-D3-Sb 

c 

subroutine vumat( 

c read only variables (unmodiflable) 

1 nblock-indirTnshrinstatevTnfieldVTnpropsilannealT 

E step„time-i total_time -idt -.cmname iCOord^mp^char^lengthT 

3 propsidensityTStrai n_inc n rel_spin_i nc -1 

M temp_old istretch^old ndef grad_old t f i eld_old -i 

5 stress_oldiState_oldTener_inter n_o ldnener_inela s_o id-. 

ti temp__new stretch_new de f grad_new-. f ield_new-. 

c write only variables (modifiable) 

7 stress_new i state_newTener_intern_newiener_inelas_new) 

c 

include ' vaba_par am . inc ' 

dimension coord_mp ( nb 1 ock i * ) 1 char_length (nblock)iprops(npro 

1 densityCnblock)-iStrain_inc(nblock-indir + nshr)n 

2 rel_sp inline Cnblock^nshr) -x temp_old ( nblock ) t 

3 stretch_old (nblockiodir+nshr) n 

M def grad_old ( nblock t ndir + nshr + nshr) ^ f ield_pld ( nblockinfieldv 

) 

S s tress_o Id (nblock ind ir + nshr) s tat e_o 1 d ( nblock instatev) 

b ener_intern_old< nblock) ^ ener_inelas_o Id (nblock ) n 

7 temp_new( nblock) s tr e tch_new C nb 1 ock t nd i r + n5hr ) -1 

a def grad_new (nblock^ndir + nshr + nshr) if i eld_new ( nblock nnfieidv 

) -» 

T s tress_new ( nblock indir + nshr) i s tate_new (nblock-«nstatev) t 
1 ener_i n t er n_new (nblock ) 1 Qner_i ne 1 as_ne w (nblock ) 

c 

integer limit 
parameter (limit=MD) 
.dimension table(2-ii) 
character*6 cmname 

parameter(zero = D.da'.one = l.dDitwo = 2.dD-ithree = 3.dDTsix = b.dQt 
1 four = ^.dO'.oneptf = l.Sd□-.zept = 0.2SdD-.twbth = D.tItltlt^btIbbbbbdDl 
E eitee=fia.dO) 

c 

c 

c props(l) - e- modulus ( temper,a ture dependent) 

c props(5) - nu- poisson ratio 

c 

c Properties 3 and M descibe the rate sensitivity of yield ba 

sed on a plot of 



c 

ps ) 
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c yield stress (x-axis) vs InCstrain rate) y-axis 

c 

c props(3) - rate sensitivity (temperature dependent) SLOPE 

c props(M) - intrinsic flow rate (temperature dependen t ) INTER 

CPT 
c 

c Property 5 descibes the pressure sensitivity of yield 

c 

c propsCS) - pressure sensitivity factor 

c 

c Property b is the fai lure cri ter ion ... either an equivalen 

t plastic strain 

c for ductile failure or a maximum principal stress for britt 

le failure 

c 

c props(ti) - failure criterion 

c 

c NOTE -THESE FOLLOWING TliiO LINES WOULD APPEAR IN THE ABACUS 

EXPLICIT 

c INPUT DECK 

c 

c *USER riATERIAL.C0NSTANTS = 5 

c 2.2MeT-.0.MD-.3.21e-7-.l.Mfie-m-.a.lb 

c *DEPVAR .DELETE=b 

c a 

c 

c 

c material properties 

emod^props ( 1 ) 
enu=props ( 5 ) 

ebulk3=emod/(one-two*enu) 
eg2=emod/(one+enu) 
sg=eg2/ two 
eg3=three*eg 
elam=(ebulk3-egE)/three 
elp2g=elam+eg2 
ratesf =props ( 3 ) 
dtebsD = dt*props (M ) 
ps f =*pr ops ( 5 ) 
.'rrates=one/ratesf 
f ailst=props (b) 
tablG(lil)=a.a 
table(E-.l)=D.D 
table(l-.E)=b.E 
table(2-.a)=0.15 
table(l'.3)=17-H3 
table(2-.3)=0.3S 
tableCliM ) =3M • M7 
table(5TH)=a.SS 
tabled. S) =53 .OT 
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table(2.S)=a.75 
tabled. b)=70.3E 

tablG(l-.7)="^i.ai 

table(2-.7)=l.lS 

table(l'.a)=14b.lti 

table<2'.a)=1.35 

table(l-.^)=aD1.3 

table(E.T)=l.SS 

c 

do lOD i = lTnblocl< 

c 

c initializQ state variables 

eqplas=state_old ( i il) 
s m_o 1 d = s t a t e_o 1 d ( i -I E ) 
icont = sta te_old ( i t3) 
tstart=total_time-dt 
if ( tstar t • 1 1 . 1 - e-ti ) then 
icont=l 

state_old ( i Tb) =one 
endi f 

c 

if (state^oldCiib) .It.D-S) then 
Stat e_n GwCi-»b)=zero 
goto 
end i f 

c 

c get hardening modulus and intrinsic resistance at t 

hard=(table(lTicont+l)-table(lTiCGnt) )/ 
1 (tablG(Eticont+l)-table(5iicont)) 

s_intr=tablG(lTicont)+hard*(eqplas-table(STicont)) 

c 

c calculate predictor stress 

traceE=5train_inc(iil) +strain_inc ( i lE ) 
del_e33=-elam*traceS/Glp2g 

s igllo = s tress_o I d ( i 1 ) + eg 5 *s t r a i n_i nc C i ^ 1) 
sig52o = s tress_old ( i . 5 ) + GgE*s tr ai n_i nc ( i -»2 ) 
sig33=zero 

sigl2=stress„old (iTM)+eg2*strain_inc(iTM) 
' . ■ ssl2s = six*( sigl2**5) 

c 

c since s t r ai n_i nc ( i 3 ) is not known apriori^ loop 3 

c times without checking for convergence (works very well 

c in practise by reducing sig33 to 0 • aDD0001*sy i e 1 d ) 

do 200 ii=li3 

tr ace = t r aceE+de 1_g3 3 

sigll=sigllo+elam*trace 

sig22 = sig22o-*'elam*tracG 

c 

c calculate equivalent von mises stress from deviatoric 
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component of trial (predictor) stress- 
smises=(sigll-sig2E)**2+(sigE5)**a+(sigll)*«E 
smises=smises+ssl5s 
smises=sqr t ( smises/ two ) 

avoid division by zero during first iteration 
smises=max(oneTSmises) 

separate the hydrostatic from the deviatoric stress 

calculate the flow direction 

shydro=(sigll+sig25) /three 

flowll= (sigll-shydro) /smises 

flowEE=(sigEE-shydro) /smises 

flow33=(sig33-shydro)/smises 

flowlE=siglE/smises 

solve for equivalent von mises stress and equivalent 

plastic strain increment 

adfp=-psf*shydro*ratesf 

deqpl=dtebsa*exp(( sm_o ld-s_intr ) *ra tesf+adfp ) 
s m_n eu = smises-eg3*deqpl 

update e33 

opfe=oneptf»deqpl 

d_epll=opf e*f 1 owll 

d_epES=opf e*f lowEE 

d_ep33=opf e*f 1 ow33 

d_eplE = op f e *f lowlE 

d_eell=s tra inline ( i 1 1) -d_epll 

d_ee5E = s tr a i n_inc ( i ^ ? ) -d_epE2 

d_ee33=-elam* ( d_eell+d_eeaE ) /elpEg 

d_eelE=strain_inc ( i T M ) -d_eplE 

del_e33=d_ee33+d_ep33 

continue 

esi = strain_inc ( i t 1) **E + s train_inc ( i i2) **5 + 
del_e33**E+two*strain_inc ( i iM) **E 
esi=sqrt(esi*twbth) 
strain_inc(iT3) =del_e33 

update stressT equivalent plastic strainn location 

of plastic strain counter and state variables 

s tre5s_new ( i n 1 ) = flowll*sm_new + s hydro 

s tres3_.new ( i i 2 ) = f 1 ow25*sm_new+shydro 

s tress_new ( i i 3 ) =zero 

s tress_new ( i t M ) = f lowlE*sm_new 

eqplas=eqplas+deqpl 

if (eqplas. gt. table (Eticont + D) icont = icont + l 
cstate_new ( i 1) = s ta te.old ( i t 1) +d_eell 
cstate_new ( i lE ) =state_old ( i lE ) +d_eeES 
cstate_new ( i 3 ) = s tate_old (113) +d_eelE 
cstate_new ( i M ) = s tate_old ( i M ) + d_epll 



/^/^- /£^ 
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cstate_new ( i 5 ) = s tate_old ( i i 5) + d_epEE 
cs ta te_new ( i t b ) =state_old (lib) +d_epl2 

save state variables: plastic strain-i vm stressi-t total 
strain ratei plastic strain ratei failure criterion flag 
5 1 a t e_neu (i'«l)=eqplas 
sta te_new ( i t E ) =sm_new 
state_new (iT3)=icont 
state_new( iTM)=esi/dt 
state_new( i-.S)=deqpl/dt 
state_new ( i -.b) = state_old ( i ib ) 

bee = - ( stress^new ( i nl) + stress_neu ( i t2 ) ) 
bee2=bee*bee 

cee = s tress_neu ( i 1 1) *stress__new ( i i 2 ) - stres5_new ( i -i M ) * 

s tr ess_new ( i -i ^ ) 

froot=bee2-four*cee 

ffrot=max(oneTfroot) 

sqbmMc=sqrt(ffrot) 

p m a X = C - b e e + s q b m M c ) / 1 w o 

p m i n = ( - b e e - s q b m ^ c ) / 1 w o 

state__new ( iT7)=pmax 

state„new( iTa)=pniin 

UNPAINTED 

failst=&S-Ob 

if (pmax.gt-failst) sta te_new ( i t b ) = zero 
strain based failure criterion 
if (eqplas.gt-failst) s tate_new (iTb)=zero 
update plastic dissipation 

p 1 as ti c_work_i nc = cleqpl * ( sm_o 1 d + sm_new ) / two 
9ner_ine las_new ( i ) =ener_inelas_oldC i ) + 
pi ast ic_work_inc/densi ty ( i ) 

continue 

return 

end 



